Introduction

Monte Carlo simulation example

In this example, we compare the power of two tests for autocorrelation:

Durbin-Watson
Breusch-Godfrey

We use two different specifications of a simulated linear regression for this purpose:

\[ y_t = \beta_0 + \beta_1 t + u_t \]

and

\[ y_t = \beta_0 + \beta_1 y_{t-1} + u_t \] where the regression coefficients of the simulated model are: \(\beta_0 = 0.25\) and \(\beta_1 = -0.75\) for both our regression specifications (deterministic trend, dynamic model). \(u_t\) is assumed stationary, with mean zero.

We want to evaluate the small-sample powers of the two tests (\(\alpha = 0.05\)) under the following scenarios:

Autocorrelation: \(\rho = 0, 0.2, 0.4, 0.6, 0.8, 0.9, 0.95, 0.99\)
Sample size: \(n = 15, 30, 50\).

Data:

For this Monte Carlo example, we need to define the DGP process first.

dgp <- function(nobs = 15, model = c("trend", "dynamic"),
                corr = 0, coef = c(0.25, -0.75), sd = 1)
{
model <- match.arg(model)
coef <- rep(coef, length.out = 2)
err <- as.vector(filter(rnorm(nobs, sd = sd), corr, method = "recursive"))
if(model == "trend") {
   x <- 1:nobs
   y <- coef[1] + coef[2] * x + err
} else {
   y <- rep(NA, nobs)
   y[1] <- coef[1] + err[1]
   for(i in 2:nobs)
     y[i] <- coef[1] + coef[2] * y[i-1] + err[i]
     x <- c(0, y[1:(nobs-1)])
}
return(data.frame(y = y, x = x))
}

Basic data plot:

(one example of simulated data, \(n=50\))

df1 <- dgp(model="trend", nobs = 50)
plot(df1$y~df1$x,main="Trend")

df2 <- dgp(model="dynamic", nobs = 50)
plot(df2$y~df2$x,main="Dynamic")

Next, we define a function that evaluates test powers:

simpower <- function(nrep = 100, size = 0.05, ...)
{
pval <- matrix(rep(NA, 2 * nrep), ncol = 2)
colnames(pval) <- c("dwtest", "bgtest")
for(i in 1:nrep) {
  dat <- dgp(...)
  pval[i,1] <- dwtest(y ~ x, data = dat,
  alternative = "two.sided")$p.value
  pval[i,2] <- bgtest(y ~ x, data = dat)$p.value
}
return(colMeans(pval < size))
}

Higher precision of this function may be obtained by increasing the nrep argument. However, computation time may be a prohibitive factor here.

We prepare a function that does the actual Monte Carlo simulation task:

simulation <- function(corr = c(0, 0.2, 0.4, 0.6, 0.8, 0.9, 0.95, 0.99), 
                       nobs = c(15, 30, 50), 
                       model = c("trend", "dynamic"), ...)
{
prs <- expand.grid(corr = corr, nobs = nobs, model = model)
nprs <- nrow(prs)
pow <- matrix(rep(NA, 2 * nprs), ncol = 2)
for(i in 1:nprs) pow[i,] <- simpower(corr = prs[i,1],
   nobs = prs[i,2], model = as.character(prs[i,3]), ...)
rval <- rbind(prs, prs)
rval$test <- factor(rep(1:2, c(nprs, nprs)),
labels = c("dwtest", "bgtest"))
rval$power <- c(pow[,1], pow[,2])
rval$nobs <- factor(rval$nobs)
return(rval)
}

Now, we can run the simulation:

set.seed(987)
psim <- simulation()
str(psim)

## 'data.frame':    96 obs. of  5 variables:
##  $ corr : num  0 0.2 0.4 0.6 0.8 0.9 0.95 0.99 0 0.2 ...
##  $ nobs : Factor w/ 3 levels "15","30","50": 1 1 1 1 1 1 1 1 2 2 ...
##  $ model: Factor w/ 2 levels "trend","dynamic": 1 1 1 1 1 1 1 1 1 1 ...
##  $ test : Factor w/ 2 levels "dwtest","bgtest": 1 1 1 1 1 1 1 1 1 1 ...
##  $ power: num  0.02 0.06 0.16 0.38 0.56 0.61 0.59 0.64 0.03 0.09 ...
##  - attr(*, "out.attrs")=List of 2
##   ..$ dim     : Named int [1:3] 8 3 2
##   .. ..- attr(*, "names")= chr [1:3] "corr" "nobs" "model"
##   ..$ dimnames:List of 3
##   .. ..$ corr : chr [1:8] "corr=0.00" "corr=0.20" "corr=0.40" "corr=0.60" ...
##   .. ..$ nobs : chr [1:3] "nobs=15" "nobs=30" "nobs=50"
##   .. ..$ model: chr [1:2] "model=trend" "model=dynamic"

Simulation output can be summarized and plotted:

Test powers, i.e. \(H_0\) rejection probabilities for increasing \(\rho\) can be observed from the table and from plots:

tab <- xtabs(power ~ corr + test + model + nobs, data = psim)
ftable(tab, row.vars = c("model", "nobs", "test"), col.vars = "corr")

##                     corr    0  0.2  0.4  0.6  0.8  0.9 0.95 0.99
## model   nobs test                                               
## trend   15   dwtest      0.02 0.06 0.16 0.38 0.56 0.61 0.59 0.64
##              bgtest      0.04 0.02 0.05 0.13 0.29 0.28 0.28 0.41
##         30   dwtest      0.03 0.09 0.46 0.77 0.92 0.99 0.96 0.97
##              bgtest      0.04 0.08 0.33 0.69 0.86 0.97 0.89 0.97
##         50   dwtest      0.06 0.24 0.72 0.95 1.00 1.00 1.00 1.00
##              bgtest      0.08 0.20 0.60 0.88 0.99 1.00 1.00 1.00
## dynamic 15   dwtest      0.00 0.01 0.00 0.01 0.01 0.01 0.01 0.01
##              bgtest      0.06 0.04 0.03 0.07 0.16 0.22 0.19 0.26
##         30   dwtest      0.02 0.00 0.00 0.00 0.03 0.01 0.04 0.13
##              bgtest      0.09 0.03 0.14 0.41 0.55 0.57 0.60 0.74
##         50   dwtest      0.01 0.01 0.03 0.01 0.05 0.14 0.31 0.69
##              bgtest      0.06 0.08 0.50 0.72 0.87 0.91 0.91 0.99

xyplot(power ~ corr | model + nobs, groups = ~ test, data = psim, type = "b")

Durbin-Watson test (blue) performs better in the trend model - the difference is only pronounced for smaller samples and lower \(\rho\).
In the dynamic model, Durbin-Watson test has almost no power except for very high correlations.
Breusch-Godfrey test (pink) performs reasonably well, unless \(n\) is very small.

Example amended from the textbook Gentleman, Hornik, Parmigiani: Use R!
For detailed discussion see chapter 7 of Use R!.